# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in sex while accounting for sex differences in maximum rufous band score.ind_prop_molt_scores <- dat %>%mutate(score =as.numeric(score)) %>%left_join(., select(ind_breeding_scores, -sex), by =c("rings_comb","season")) %>%filter(!is.na(max_breeding_score)) %>%mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%select(rings_comb, season, date_J, sex, score, max_breeding_score, prop_molt_score) %>%arrange(rings_comb, season, date_J) %>%distinct() %>%filter(!is.na(prop_molt_score))# Assess sample sizes of each sexind_prop_molt_scores %>%group_by(sex) %>%summarise(n_distinct(rings_comb))
# A tibble: 2 × 2
sex `n_distinct(rings_comb)`
<chr> <int>
1 F 49
2 M 35
Code
# Step 1 and 2: Standardize max_breeding_score within each sex groupind_prop_molt_scores <- ind_prop_molt_scores %>%group_by(sex) %>%mutate(std_max_breeding_score = (max_breeding_score -mean(max_breeding_score, na.rm =TRUE)) /sd(max_breeding_score, na.rm =TRUE)) %>%ungroup()# mixed effects binomial model comparing sex and date effect on the changes in moult scoresmod1 <- lme4::glmer(prop_molt_score ~ date_J * sex + std_max_breeding_score +# date_J * sex + (1| rings_comb) + (1|season),data = ind_prop_molt_scores, family = binomial,control =glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =1e+06)))plot(allEffects(mod1))
Code
# strong date effect, but no sex effecttbl_regression(mod1, intercept =TRUE, label =list(date_J ~"Date", sex ~"Sex"))
Characteristic
log(OR)1
95% CI1
p-value
(Intercept)
20
15, 24
<0.001
Date
-0.09
-0.11, -0.07
<0.001
Sex
F
—
—
M
1.1
-5.7, 8.0
0.8
std_max_breeding_score
-0.29
-0.65, 0.07
0.11
Date * Sex
Date * M
-0.01
-0.04, 0.03
0.7
1 OR = Odds Ratio, CI = Confidence Interval
Code
# extract predicted trendspr <- ggeffects::predict_response(mod1, c("date_J [30:293]", "sex"))# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian datedates_for_plot <-data.frame(date =as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),date_J =c(1:365))# join the back-transformed dates to model fitsmod1_fits <-as.data.frame(pr) %>%rename(date_J = x,sex = group) %>%left_join(., dates_for_plot, by ="date_J")# plot the modelggplot() +geom_line(data = mod1_fits, aes(x = date, y = predicted, color = sex)) +geom_ribbon(data = mod1_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = sex),lwd =1, alpha =0.25) + luke_theme +theme(legend.position =c(0.3, 0.2),legend.justification =c(1, 0),strip.background =element_blank(),axis.title.x =element_blank(),axis.text.x =element_text(size =10, angle =45, hjust =1, vjust =1)) +ylab("proportion intact of individual-based maximum rufous band") +scale_x_date(date_labels ="%B", expand =c(0.01, 0.01), date_breaks ="1 month") +scale_colour_brewer(palette ="Dark2", direction =-1,name ="Sex",labels =c("Female (N = 49)", "Male (N = 35)")) +scale_fill_brewer(palette ="Dark2", direction =-1,name ="Sex",labels =c("Female (N = 49)", "Male (N = 35)"))
# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in migratory status.ind_prop_molt_scores <- dat %>%mutate(score =as.numeric(score)) %>%left_join(., select(ind_breeding_scores, -sex), by =c("rings_comb","season")) %>%filter(!is.na(max_breeding_score)) %>%mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%select(rings_comb, season, date_J, sex, migratory_status, score, max_breeding_score, prop_molt_score) %>%arrange(rings_comb, season, date_J) %>%distinct() %>%filter(!is.na(prop_molt_score)) %>%filter(migratory_status %in%c("R", "M"))# Assess sample sizes of each sexind_prop_molt_scores %>%group_by(sex, migratory_status) %>%summarise(n_distinct(rings_comb))
# A tibble: 4 × 3
# Groups: sex [2]
sex migratory_status `n_distinct(rings_comb)`
<chr> <chr> <int>
1 F M 2
2 F R 11
3 M M 3
4 M R 5
Code
# mixed effects binomial model comparing sex and date effect on the changes in moult scoresmod1_status <- lme4::glmer(prop_molt_score ~ date_J * migratory_status +# date_J * sex + (1| rings_comb) + (1| season),data = ind_prop_molt_scores, family = binomial,control =glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =1e+06)))# strong date effect, but no status effecttbl_regression(mod1_status, intercept =TRUE, label =list(date_J ~"Date", migratory_status ~"Migratory status"))
Characteristic
log(OR)1
95% CI1
p-value
(Intercept)
25
10, 40
0.001
Date
-0.11
-0.19, -0.04
0.002
Migratory status
M
—
—
R
12
-17, 41
0.4
Date * Migratory status
Date * R
-0.06
-0.20, 0.08
0.4
1 OR = Odds Ratio, CI = Confidence Interval
Code
# extract predicted trendspr <- ggeffects::predict_response(mod1_status, c("date_J [30:293]", "migratory_status"))# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian datedates_for_plot <-data.frame(date =as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),date_J =c(1:365))# join the back-transformed dates to model fitsmod1_fits <-as.data.frame(pr) %>%rename(date_J = x,migratory_status = group) %>%left_join(., dates_for_plot, by ="date_J")# plot the modelggplot() +geom_line(data = mod1_fits, aes(x = date, y = predicted, color = migratory_status)) +geom_ribbon(data = mod1_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = migratory_status),lwd =1, alpha =0.25) + luke_theme +theme(legend.position =c(0.3, 0.2),legend.justification =c(1, 0),strip.background =element_blank(),axis.title.x =element_blank(),axis.text.x =element_text(size =10, angle =45, hjust =1, vjust =1)) +ylab("proportion intact of individual-based maximum rufous band") +scale_x_date(date_labels ="%B", expand =c(0.01, 0.01), date_breaks ="1 month") +scale_colour_brewer(palette ="Set1", direction =-1,name ="Migratory status",labels =c("Resident (N = 16)", "Migratory (N = 5)")) +scale_fill_brewer(palette ="Set1", direction =-1,name ="Migratory status",labels =c("Resident (N = 16)", "Migratory (N = 5)"))
within-individual moult dynamics (age)
Code
# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in migratory status.ind_prop_molt_scores <- dat %>%mutate(score =as.numeric(score)) %>%left_join(., select(ind_breeding_scores, -sex), by =c("rings_comb","season")) %>%filter(!is.na(max_breeding_score)) %>%mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%select(rings_comb, season, date_J, sex, age, score, max_breeding_score, prop_molt_score) %>%arrange(rings_comb, season, date_J) %>%distinct() %>%filter(!is.na(prop_molt_score)) %>%filter(!is.na(age))# Assess sample sizes of each sexind_prop_molt_scores %>%group_by(age) %>%summarise(n_distinct(rings_comb))
# mixed effects binomial model comparing sex and date effect on the changes in moult scoresmod1_age <- lme4::glmer(prop_molt_score ~ date_J + age +# date_J * sex + (1| rings_comb),# + (1 | season),data = ind_prop_molt_scores, family = binomial)# control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))# strong date effect, but no status effecttbl_regression(mod1_age, intercept =TRUE, label =list(date_J ~"Date", age ~"Age"))
Characteristic
log(OR)1
95% CI1
p-value
(Intercept)
53
18, 88
0.003
Date
-0.28
-0.47, -0.10
0.003
Age
2.2
0.61, 3.7
0.006
1 OR = Odds Ratio, CI = Confidence Interval
Code
# extract predicted trendspr <- ggeffects::predict_response(mod1_age, c("date_J [30:293]", "age"))# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian datedates_for_plot <-data.frame(date =as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),date_J =c(1:365))# join the back-transformed dates to model fitsmod1_fits <-as.data.frame(pr) %>%rename(date_J = x,age = group) %>%left_join(., dates_for_plot, by ="date_J")# plot the modelggplot() +geom_line(data = mod1_fits, aes(x = date, y = predicted, color = age)) +geom_ribbon(data = mod1_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = age),lwd =1, alpha =0.25) + luke_theme +theme(legend.position =c(0.3, 0.2),legend.justification =c(1, 0),strip.background =element_blank(),axis.title.x =element_blank(),axis.text.x =element_text(size =10, angle =45, hjust =1, vjust =1)) +ylab("proportion intact of individual-based maximum rufous band") +scale_x_date(date_labels ="%B", expand =c(0.01, 0.01), date_breaks ="1 month", limits =c(as.Date("2021-10-01"), as.Date("2022-05-01"))) +scale_colour_brewer(palette ="Spectral", direction =-1,name ="Age",labels =c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", "4 (N = 13", "5 (N = 4)", "6 (N = 2)")) +scale_fill_brewer(palette ="Spectral", direction =-1,name ="Age",labels =c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", "4 (N = 13", "5 (N = 4)", "6 (N = 2)"))
13-Nov-24: meeting notes (Bart, Luke, and Bashar via zoom)
add the breeding data to the analysis: for each individual include the date of last breeding (i.e., last date seen with nest or brood) as a predictor of molt timing. Also do an analysis relating the timing of last breeding to the age of the individual. Bashar will send Luke the breeding data he received from Ted.
drop pre-October data from the molt timing analysis.
next steps: Bashar to eventually write up this study as a report (which could eventually be used as his MSc thesis if he wants). Priority will be on him writing up a proposal for his thesis (which has a concrete deadline and is graded). Bashar will present the study at our weekly Monday meeting in the first few months of 2025.
Source Code
---title: "Banded Dotterel Moult Study"subtitle: | Exploration of datasetdate: "`r format(Sys.time(), '%d %B, %Y')`"author: - name: Luke Eberhart-Hertel orcid: 0000-0001-7311-6088 email: luke.eberhart@bi.mpg.de url: https://www.bi.mpg.de/person/115852/2867 affiliations: - ref: bk - name: Bashar Jarayseh affiliations: - ref: bk - name: Ailsa Howard affiliations: - ref: ah - name: Tony Habraken affiliations: - ref: th - name: Emma Williams affiliations: - ref: ew - name: Colin O`Donnell affiliations: - ref: ew - name: Bart Kempenaers affiliations: - ref: bkaffiliations: - id: bk number: 1 name: Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany - id: ah number: 2 name: South Bay Banded Dotterel Project, Kaikoura, New Zealand - id: th number: 2 name: Port Waikato Banded Dotterel Project, Port Waikato, New Zealand - id: ew number: 3 name: Department of Conservation, Christchurch, New Zealandformat: html: toc: true code-fold: true code-tools: true self-contained: true highlight-style: github theme: Cosmoexecute: warning: false cache: trueeditor_options: chunk_output_type: console---```{r}knitr::opts_chunk$set(cache =TRUE)``````{r, message=FALSE, results='hide', warning=FALSE, results='hide', cache=FALSE, include=FALSE}## Prerequisites### R packages# - The following packages are needed for analysis and can be easily installed from [CRAN](http://cran.r-project.org/) or GitHub by running the following code chunk:# a vector of all the packages needed in the projectpackages_required_in_project <-c("tidyverse","readxl","RMark","RColorBrewer","patchwork","mapview","lubridate","extrafont","here","DT","leaflet","sf","leafpop","tsibble","corrplot","gghalves","gam","pscl","gamlss","gt","lme4","ggpattern","gtsummary","effects","lattice","rptR","partR2","broom.mixed","forcats","glmmTMB")# of the required packages, check if some need to be installednew.packages <- packages_required_in_project[!(packages_required_in_project %in%installed.packages()[,"Package"])]# install all packages that are not locally availableif(length(new.packages)) install.packages(new.packages)# load all the packages into the current R sessionlapply(packages_required_in_project, require, character.only =TRUE)# set the home directory to where the project is locally based (i.e., to find # the relevant datasets to import, etc.here::set_here()``````{r, message=FALSE, results='hide', warning=FALSE, include=FALSE}### Plotting themes# - The following plotting themes, colors, and typefaces are used throughout the project:# Find fonts from computer that you want. Use regular expressions to do this# For example, load all fonts that are 'verdana' or 'Verdana'extrafont::font_import(pattern ="[V/v]erdana", prompt =FALSE) # check which fonts were loadedextrafont::fonts()extrafont::fonttable()extrafont::loadfonts() # load these into R# define the plotting theme to be used in subsequent ggplotsluke_theme <-theme_bw() +theme(text =element_text(family ="Verdana"),legend.title =element_text(size =10),legend.text =element_text(size =8),axis.title.x =element_text(size =10),axis.text.x =element_text(size =8), axis.title.y =element_text(size =10),axis.text.y =element_text(size =8),strip.text =element_text(size =10),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),axis.ticks =element_line(linewidth =0.5, colour ="grey40"),axis.ticks.length =unit(0.2, "cm"),panel.border =element_rect(linetype ="solid", colour ="grey"),legend.position.inside =c(0.1, 0.9) )```## Explore Bashar's Dataset (modified 7-Nov-2024)This dataset contains the scored moult data for all usable photos from the 2021-2022, 2022-2023, and 2023-2024 breeding seasons at Kaikoura```{r, include = FALSE}# import data with all columns as character (so that no auto-formatting is done by R)dat <-read.csv(here("../data/Dataset_Molt_Final - modified 7.11.2024.csv"), colClasses ="character") %>%mutate(Rings_comb =str_sub(Rings_comb, 1, 4)) %>%rename_with(tolower) # convert all column headers to lowercase# check values of Date column for mistakes...looks goodunique(dat$date)# check values of score column for mistakes...one row to fixunique(dat$score)# one row contains no data for the score (row 767)...needs to be fixeddat %>%filter(score =="")# check values of rings_comb column for mistakes...need to remove the non-unique combos (rows with XB and XW)unique(dat$rings_comb)dat %>%filter(rings_comb %in%c("XB", "XW"))dat <- dat %>%filter(!grepl("XB", rings_comb) &!grepl("XW", rings_comb))# import sex, migratory status, and age at banding informationdat_sexes <-read.csv(here("../data/sex_status_banding.csv"), colClasses ="character") %>%mutate(Rings_comb =str_sub(Rings_comb, 1, 4)) %>%rename_with(tolower) %>%# convert all column headers to lowercasedistinct()# checkout the sex-type datadat_sexes %>%pull(rings_comb) %>%unique() %>%length()dat_sexes[which(duplicated(dat_sexes)), ]dat_sexes %>%filter(migratory_status %in%c("R", "M"))dat_sexes %>%filter(rings_comb =="RBYG")# mutate the Date column into a date variabledat <- dat %>%mutate(date =paste(substring(dat$date, first =7, last =10), substring(date, first =4, last =5),substring(date, first =1, last =2),sep ="-") %>%as.Date()) %>%# subset to data with Molt == 1filter(molt ==1) %>%# specify the season as the first calender yearmutate(season =ifelse(month(date) <7, year(date) -1, year(date))) %>%# change to Julian date shifted for the Southern Hemisphere (1 = July 1)mutate(date_J =as.numeric(format(date +181, "%j"))) %>%# join the sexes provided by Ailsa (note many-to-many is fine)left_join(., dat_sexes, by ="rings_comb", relationship ="many-to-many") %>%# remove individuals with unknown sexfilter(sex !="U") %>%filter(sex !="N") %>%# add a ranking variable to sort the facets of the sampling distribution plotsgroup_by(rings_comb) %>%mutate(n_photos =n(), n_scores =n_distinct(score)) %>%arrange(desc(n_scores)) %>%mutate(banding_date =paste(substring(banding_date, first =7, last =10), substring(banding_date, first =4, last =5),substring(banding_date, first =1, last =2),sep ="-") %>%as.Date(),season_ringed =ifelse(age_at_banding =="P" , ifelse(month(banding_date) <7, year(banding_date) -1, year(banding_date)),NA),age =ifelse(age_at_banding =="P", as.numeric(season - season_ringed), NA))# i should find a better way to calculate agedat %>%group_by(season) %>%summarise(n_obs =n(),min_date =min(date),max_date =max(date))# summarize the number of individuals in the dataset...95 unique combosdat %>%pull(rings_comb) %>%unique() %>%length()# summarise the number of seasons for each individual...63 have more than one season of datadat %>%select(rings_comb, season) %>%distinct() %>%group_by(rings_comb) %>%summarise(n_seasons =n()) %>%arrange(desc(n_seasons)) %>%filter(n_seasons >1) %>%nrow()# summarise the number of usable photos for each individual...dat %>%group_by(rings_comb, season) %>%summarise(n_photos =n()) %>%arrange(desc(n_photos))```### Plot the sampling distribution of Ailsa's photos across the seasons```{r, fig.height=9}# 2021-2022 season dataggplot(data = dat %>%filter(season ==2021& sex =="F") %>%mutate(score =as.numeric(score), rings_comb_ =reorder(rings_comb, n_photos, decreasing =TRUE))) +geom_point(aes(y =1, x = date, fill = score), pch =21, color ="black", size =3) +facet_wrap(. ~ rings_comb_, ncol =1, strip.position ="right") +scale_fill_gradient(high ="#cc4c02", low ="white") +theme_bw() +theme(axis.title.y =element_blank(),axis.text.y =element_blank(),panel.grid.major.y =element_blank(),panel.grid.minor =element_blank(),axis.ticks =element_blank(),strip.text.y.right =element_text(angle =0)) +scale_x_date(date_labels ="%W", expand =c(0.01, 0.01), date_breaks ="3 week", limits =c(as.Date("2021-07-05"), as.Date("2022-05-01"))) +xlab("week") +ggtitle("Females from the 2021-2022 breeding season in Kaikoura")``````{r, fig.height=8}ggplot(data = dat %>%filter(season ==2021& sex =="M") %>%mutate(score =as.numeric(score))) +geom_point(aes(y =1, x = date, fill = score), pch =21, color ="black", size =3) +facet_wrap(. ~reorder(rings_comb, n_photos, decreasing =TRUE), ncol =1, strip.position ="right") +scale_fill_gradient(high ="#cc4c02", low ="white") +theme_bw() +theme(axis.title.y =element_blank(),axis.text.y =element_blank(),panel.grid.major.y =element_blank(),panel.grid.minor =element_blank(),axis.ticks =element_blank(),strip.text.y.right =element_text(angle =0)) +scale_x_date(date_labels ="%W", expand =c(0.01, 0.01), date_breaks ="3 week", limits =c(as.Date("2021-07-05"), as.Date("2022-05-01"))) +xlab("week") +ggtitle("Males from the 2021-2022 breeding season in Kaikoura")``````{r, fig.height=10}# 2022-2023 season data (not as good coverage as the previous season)ggplot(data = dat %>%filter(season ==2022& sex =="F") %>%mutate(score =as.numeric(score))) +geom_point(aes(y =1, x = date, fill = score), pch =21, color ="black", size =3) +facet_wrap(. ~reorder(rings_comb, n_photos, decreasing =TRUE), ncol =1, strip.position ="right") +scale_fill_gradient(high ="#cc4c02", low ="white") +theme_bw() +theme(axis.title.y =element_blank(),axis.text.y =element_blank(),panel.grid.major.y =element_blank(),panel.grid.minor =element_blank(),axis.ticks =element_blank(),strip.text.y.right =element_text(angle =0)) +scale_x_date(date_labels ="%W", expand =c(0.01, 0.01), date_breaks ="3 week", limits =c(as.Date("2022-07-05"), as.Date("2023-05-01"))) +xlab("week") +ggtitle("Females from the 2022-2023 breeding season in Kaikoura")``````{r, fig.height=7}ggplot(data = dat %>%filter(season ==2022& sex =="M") %>%mutate(score =as.numeric(score))) +geom_point(aes(y =1, x = date, fill = score), pch =21, color ="black", size =3) +facet_wrap(. ~reorder(rings_comb, n_photos, decreasing =TRUE), ncol =1, strip.position ="right") +scale_fill_gradient(high ="#cc4c02", low ="white") +theme_bw() +theme(axis.title.y =element_blank(),axis.text.y =element_blank(),panel.grid.major.y =element_blank(),panel.grid.minor =element_blank(),axis.ticks =element_blank(),strip.text.y.right =element_text(angle =0)) +scale_x_date(date_labels ="%W", expand =c(0.01, 0.01), date_breaks ="3 week", limits =c(as.Date("2022-07-05"), as.Date("2023-05-01"))) +xlab("week") +ggtitle("Males from the 2022-2023 breeding season in Kaikoura")``````{r, fig.height=10}# 2023-2024 season data ggplot(data = dat %>%filter(season ==2023& sex =="F") %>%mutate(score =as.numeric(score))) +geom_point(aes(y =1, x = date, fill = score), pch =21, color ="black", size =3) +facet_wrap(. ~reorder(rings_comb, n_photos, decreasing =TRUE), ncol =1, strip.position ="right") +scale_fill_gradient(high ="#cc4c02", low ="white") +theme_bw() +theme(axis.title.y =element_blank(),axis.text.y =element_blank(),panel.grid.major.y =element_blank(),panel.grid.minor =element_blank(),axis.ticks =element_blank(),strip.text.y.right =element_text(angle =0)) +scale_x_date(date_labels ="%W", expand =c(0.01, 0.01), date_breaks ="3 week", limits =c(as.Date("2023-07-01"), as.Date("2024-05-01"))) +xlab("week") +ggtitle("Females from the 2023-2024 breeding season in Kaikoura")``````{r, fig.height=6}ggplot(data = dat %>%filter(season ==2023& sex =="M") %>%mutate(score =as.numeric(score))) +geom_point(aes(y =1, x = date, fill = score), pch =21, color ="black", size =3) +facet_wrap(. ~reorder(rings_comb, n_photos, decreasing =TRUE), ncol =1, strip.position ="right") +scale_fill_gradient(high ="#cc4c02", low ="white") +theme_bw() +theme(axis.title.y =element_blank(),axis.text.y =element_blank(),panel.grid.major.y =element_blank(),panel.grid.minor =element_blank(),axis.ticks =element_blank(),strip.text.y.right =element_text(angle =0)) +scale_x_date(date_labels ="%W", expand =c(0.01, 0.01), date_breaks ="3 week", limits =c(as.Date("2023-07-05"), as.Date("2024-05-01"))) +xlab("week") +ggtitle("Males from the 2023-2024 breeding season in Kaikoura")```### sex-specific differences in breeding plumage```{r, include=FALSE}# extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual.ind_breeding_scores <- dat %>%mutate(score =as.numeric(score)) %>%mutate(breeding_season =ifelse(month(date) %in%c(8, 9, 10, 11, 12), 1, 0)) %>%filter(breeding_season ==1) %>%group_by(season, rings_comb, sex) %>%summarise(max_breeding_score =max(score)) %>%ungroup()# Assess sample sizes of each sex (49 females, 35 males)ind_breeding_scores %>%group_by(sex) %>%summarise(n_distinct(rings_comb))# linear mixed model for the difference in max score between sexesmod_max_score <-lmer(max_breeding_score ~ sex + (1| rings_comb) + (1| season), data = ind_breeding_scores)summary(mod_max_score)plot(allEffects(mod_max_score))# Derive confidence intervals of effect sizes from parametric bootstrappingtidy_mod_max_score <-tidy(mod_max_score, conf.int =TRUE, conf.method ="boot", nsim =1000)tbl_regression(mod_max_score, intercept =TRUE, label =list(sex ~"Sex"))# run rptR to obtain repeatabilities of random effectsrpt_mod_max_score <-rpt(max_breeding_score ~ sex + (1| rings_comb) + (1| season),grname =c("rings_comb", "season", "Fixed"),data = ind_breeding_scores,datatype ="Gaussian",nboot =1000, npermut =1000, ratio =TRUE,adjusted =TRUE, ncores =4, parallel =TRUE)# run partR2 on each model to obtain marginal R2, parameter estimates, and beta# weightsR2m_mod_max_score <-partR2(mod_max_score,partvars =c("sex"),R2_type ="marginal",nboot =1000,CI =0.95,max_level =1)R2c_mod_max_score <-partR2(mod_max_score,partvars =c("sex"),R2_type ="conditional",nboot =1000,CI =0.95,max_level =1)stats_mod_max_score <-list(mod = mod_max_score,tidy = tidy_mod_max_score,rptR = rpt_mod_max_score,partR2m = R2m_mod_max_score,partR2c = R2c_mod_max_score,data = ind_breeding_scores)# saveRDS(stats_mod_max_score, file = "out/stats_mod_max_score.rds")#### Table of effect sizes ----# Retrieve sample sizessample_sizes <- ind_breeding_scores %>%ungroup() %>%summarise(Year =n_distinct(season),Individual =n_distinct(rings_comb),Observations =nrow(.))sample_sizes <-as.data.frame(t(as.data.frame(sample_sizes))) %>%rownames_to_column("term") %>%rename(estimate = V1) %>%mutate(stat ="n")# clean model component namesmod_comp_names <-data.frame(comp_name =c("Male (Sex)","Total Marginal \U1D479\U00B2","Sex","Total Conditional \U1D479\U00B2","Individual","Year","Residual","Individual","Year","Residual","Years","Individuals","Observations"))# Fixed effect sizes (non-standardized)fixefTable <- stats_mod_max_score$tidy %>% dplyr::filter(effect =="fixed") %>% dplyr::select(term, estimate, conf.low, conf.high) %>%as.data.frame() %>%mutate(stat ="fixed")# Fixed effect sizes (standardized)fixef_bw_Table <- stats_mod_max_score$partR2m$BW %>%as.data.frame() %>%mutate(stat ="fixed_bw") %>%rename(conf.low = CI_lower,conf.high = CI_upper)# Semi-partial R2 estimatesR2Table <-bind_rows(stats_mod_max_score$partR2m$R2, stats_mod_max_score$partR2c$R2[1,]) %>% dplyr::select(term, estimate, CI_lower, CI_upper) %>%as.data.frame() %>%mutate(stat ="partR2") %>%rename(conf.low = CI_lower,conf.high = CI_upper)# Random effects variancesranefTable <- stats_mod_max_score$tidy %>% dplyr::filter(effect =="ran_pars") %>% dplyr::select(group, estimate, conf.low, conf.high) %>%as.data.frame() %>%mutate(stat ="rand") %>%rename(term = group) %>%mutate(estimate = estimate^2,conf.high = conf.high^2,conf.low = conf.low^2)# Adjusted repeatabilitiescoefRptTable <- stats_mod_max_score$rptR$R_boot %>% dplyr::select(-Fixed) %>%mutate(residual =1-rowSums(.)) %>%apply(., 2, function(x) c(mean (x), quantile (x, prob =c(0.025, 0.975)))) %>%t() %>%as.data.frame() %>%rownames_to_column("term") %>%rename(estimate = V1,conf.low =`2.5%`,conf.high =`97.5%`) %>%mutate(stat ="RptR")# Store all parameters into a single table and clean it upallCoefs_mod <-bind_rows(fixef_bw_Table, R2Table, ranefTable, coefRptTable, sample_sizes) %>%bind_cols(., mod_comp_names) %>%mutate(coefString =ifelse(!is.na(conf.low),paste0("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"),NA),effect =c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)),rep("Partitioned \U1D479\U00B2", nrow(R2Table)),rep("Random effects \U1D70E\U00B2", nrow(ranefTable)),rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)),rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>% dplyr::select(effect, everything())``````{r}# density plot of sex-specific breeding plumage score distributionsind_breeding_scores %>%group_by(rings_comb, sex) %>%summarise(max_breeding_score =max(max_breeding_score)) %>%ungroup() %>%group_by(sex, max_breeding_score) %>%summarise(n =n()) %>%bind_rows(., data.frame(sex =c("F", "M"), max_breeding_score =c(7, 3), n =c(0, 0))) %>%ungroup() %>%ggplot() +geom_bar(aes(x = max_breeding_score, y = n, fill = sex), alpha =0.5, stat ="identity", position ="dodge", width =0.5) + luke_theme +theme(legend.position =c(0.15, 0.9)) +xlab("Maximum individual rufous band score during breeding season") +ylab("Frequency") +scale_colour_brewer(palette ="Dark2", direction =-1,name ="Sex",labels =c("Female (N = 49)", "Male (N = 35)")) +scale_fill_brewer(palette ="Dark2", direction =-1,name ="Sex",labels =c("Female (N = 49)", "Male (N = 35)"))# draw gt tablemod_max_score_table <- allCoefs_mod %>% dplyr::select(effect, comp_name, estimate, coefString) %>%gt(rowname_col ="row",groupname_col ="effect") %>%cols_label(comp_name =html("<i>Banded Dotterel rufous band score</i>"),estimate ="Mean estimate",coefString ="95% confidence interval") %>%fmt_number(columns =c(estimate),rows =1:10,decimals =2,use_seps =FALSE) %>%fmt_number(columns =c(estimate),rows =11:13,decimals =0,use_seps =FALSE) %>%sub_missing(columns =1:4,missing_text ="") %>%cols_align(align ="left",columns =c(comp_name)) %>%tab_options(row_group.font.weight ="bold",row_group.background.color =brewer.pal(9,"Greys")[3],table.font.size =12,data_row.padding =3,row_group.padding =4,summary_row.padding =2,column_labels.font.size =14,row_group.font.size =12,table.width =pct(80))mod_max_score_tableplot(allEffects(mod_max_score))```### age- and sex-specific variation in breeding plumage (of known-aged birds)```{r, include=FALSE}# extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual.ind_breeding_scores_age <- dat %>%mutate(score =as.numeric(score)) %>%mutate(breeding_season =ifelse(month(date) %in%c(8, 9, 10, 11, 12), 1, 0)) %>%filter(breeding_season ==1) %>%group_by(season, rings_comb, sex, age) %>%summarise(max_breeding_score =max(score)) %>%ungroup() %>%filter(!is.na(age))# Assess sample sizes of each sex and season (14 females, 10 males)ind_breeding_scores_age %>%group_by(sex) %>%summarise(n_distinct(rings_comb))# linear mixed model for the difference in max score between sexesmod_max_score_age_sex <-lmer(max_breeding_score ~ sex + age + (1| rings_comb) + (1| season), data = ind_breeding_scores_age)summary(mod_max_score_age_sex)plot(allEffects(mod_max_score_age_sex))# Derive confidence intervals of effect sizes from parametric bootstrappingtidy_mod_max_score_age_sex <-tidy(mod_max_score_age_sex, conf.int =TRUE, conf.method ="boot", nsim =1000)tbl_regression(mod_max_score_age_sex, intercept =TRUE, label =list(sex ~"Sex", age ~"Age"))# run rptR to obtain repeatabilities of random effectsrpt_mod_max_score_age_sex <-rpt(max_breeding_score ~ sex + age + (1| rings_comb) + (1| season),grname =c("rings_comb", "season", "Fixed"),data = ind_breeding_scores_age,datatype ="Gaussian",nboot =1000, npermut =1000, ratio =TRUE,adjusted =TRUE, ncores =4, parallel =TRUE)# run partR2 on each model to obtain marginal R2, parameter estimates, and beta# weightsR2m_mod_max_score_age_sex <-partR2(mod_max_score_age_sex,partvars =c("sex", "age"),R2_type ="marginal",nboot =1000,CI =0.95,max_level =1)R2c_mod_max_score_age_sex <-partR2(mod_max_score_age_sex,partvars =c("sex", "age"),R2_type ="conditional",nboot =1000,CI =0.95,max_level =1)stats_mod_max_score_age_sex <-list(mod = mod_max_score_age_sex,tidy = tidy_mod_max_score_age_sex,rptR = rpt_mod_max_score_age_sex,partR2m = R2m_mod_max_score_age_sex,partR2c = R2c_mod_max_score_age_sex,data = ind_breeding_scores_age)# saveRDS(stats_mod_max_score_age_sex, file = "out/stats_mod_max_score_age_sex.rds")#### Table of effect sizes ----# Retrieve sample sizessample_sizes <- ind_breeding_scores_age %>%ungroup() %>%summarise(Year =n_distinct(season),Individual =n_distinct(rings_comb),Observations =nrow(.))sample_sizes <-as.data.frame(t(as.data.frame(sample_sizes))) %>%rownames_to_column("term") %>%rename(estimate = V1) %>%mutate(stat ="n")# clean model component namesmod_comp_names <-data.frame(comp_name =c("Male (Sex)","Age","Total Marginal \U1D479\U00B2","Sex","Age","Total Conditional \U1D479\U00B2","Individual","Year","Residual","Individual","Year","Residual","Years","Individuals","Observations"))# Fixed effect sizes (non-standardized)fixefTable <- stats_mod_max_score_age_sex$tidy %>% dplyr::filter(effect =="fixed") %>% dplyr::select(term, estimate, conf.low, conf.high) %>%as.data.frame() %>%mutate(stat ="fixed")# Fixed effect sizes (standardized)fixef_bw_Table <- stats_mod_max_score_age_sex$partR2m$BW %>%as.data.frame() %>%mutate(stat ="fixed_bw") %>%rename(conf.low = CI_lower,conf.high = CI_upper)# Semi-partial R2 estimatesR2Table <-bind_rows(stats_mod_max_score_age_sex$partR2m$R2, stats_mod_max_score_age_sex$partR2c$R2[1,]) %>% dplyr::select(term, estimate, CI_lower, CI_upper) %>%as.data.frame() %>%mutate(stat ="partR2") %>%rename(conf.low = CI_lower,conf.high = CI_upper)# Random effects variancesranefTable <- stats_mod_max_score_age_sex$tidy %>% dplyr::filter(effect =="ran_pars") %>% dplyr::select(group, estimate, conf.low, conf.high) %>%as.data.frame() %>%mutate(stat ="rand") %>%rename(term = group) %>%mutate(estimate = estimate^2,conf.high = conf.high^2,conf.low = conf.low^2)# Adjusted repeatabilitiescoefRptTable <- stats_mod_max_score_age_sex$rptR$R_boot %>% dplyr::select(-Fixed) %>%mutate(residual =1-rowSums(.)) %>%apply(., 2, function(x) c(mean (x), quantile (x, prob =c(0.025, 0.975)))) %>%t() %>%as.data.frame() %>%rownames_to_column("term") %>%rename(estimate = V1,conf.low =`2.5%`,conf.high =`97.5%`) %>%mutate(stat ="RptR")# Store all parameters into a single table and clean it upallCoefs_mod <-bind_rows(fixef_bw_Table, R2Table, ranefTable, coefRptTable, sample_sizes) %>%bind_cols(., mod_comp_names) %>%mutate(coefString =ifelse(!is.na(conf.low),paste0("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"),NA),effect =c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)),rep("Partitioned \U1D479\U00B2", nrow(R2Table)),rep("Random effects \U1D70E\U00B2", nrow(ranefTable)),rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)),rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>% dplyr::select(effect, everything())``````{r}# draw gt tablemod_max_score_age_sex_table <- allCoefs_mod %>% dplyr::select(effect, comp_name, estimate, coefString) %>%gt(rowname_col ="row",groupname_col ="effect") %>%cols_label(comp_name =html("<i>Banded Dotterel rufous band score</i>"),estimate ="Mean estimate",coefString ="95% confidence interval") %>%fmt_number(columns =c(estimate),rows =1:12,decimals =2,use_seps =FALSE) %>%fmt_number(columns =c(estimate),rows =11:15,decimals =0,use_seps =FALSE) %>%sub_missing(columns =1:4,missing_text ="") %>%cols_align(align ="left",columns =c(comp_name)) %>%tab_options(row_group.font.weight ="bold",row_group.background.color =brewer.pal(9,"Greys")[3],table.font.size =12,data_row.padding =3,row_group.padding =4,summary_row.padding =2,column_labels.font.size =14,row_group.font.size =12,table.width =pct(80))mod_max_score_age_sex_tableplot(allEffects(mod_max_score_age_sex))```### migratory status- and sex-specific variation in breeding plumage (of known-status birds)```{r, include=FALSE}# extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual.ind_breeding_scores_status <- dat %>%mutate(score =as.numeric(score)) %>%mutate(breeding_season =ifelse(month(date) %in%c(8, 9, 10, 11, 12), 1, 0)) %>%filter(breeding_season ==1) %>%group_by(season, rings_comb, sex, age, migratory_status) %>%summarise(max_breeding_score =max(score)) %>%ungroup() %>%filter(migratory_status %in%c("M", "R"))# Assess sample sizes of each sex and season # (13 females (2 migrant, 11 resident), 8 males (3 migrant, 5 resident))ind_breeding_scores_status %>%group_by(sex, migratory_status) %>%summarise(n_distinct(rings_comb))# linear mixed model for the difference in max score between sexesmod_max_score_status_sex <-lmer(max_breeding_score ~ sex + migratory_status + (1| rings_comb), data = ind_breeding_scores_status)# # detected singularity issue in a model with (1 | season), summary showed that variance component of this random effect was zero. # # Model with one random effectmod1 <-lmer(max_breeding_score ~ sex + migratory_status + season + (1| rings_comb), data = ind_breeding_scores_status)mod2 <-lmer(max_breeding_score ~ sex + migratory_status + (1| season), data = ind_breeding_scores_status)mod3 <-lmer(max_breeding_score ~ sex + migratory_status + (1| rings_comb), data = ind_breeding_scores_status)# Compare with the full modelmod_full <-lmer(max_breeding_score ~ sex + migratory_status + (1| rings_comb) + (1| season), data = ind_breeding_scores_status)# Check the summariessummary(mod1)summary(mod2)summary(mod3)summary(mod_full)# Compare models using AICAIC(mod1, mod2, mod3, mod_full)# conclude that the term (1 | season) is not needed for statistical reasons, but also since we are not interested in population-level annual variation in the max-breeding plumage scores we can drop thissummary(mod_max_score_status_sex)plot(allEffects(mod_max_score_status_sex))# Derive confidence intervals of effect sizes from parametric bootstrappingtidy_mod_max_score_status_sex <-tidy(mod_max_score_status_sex, conf.int =TRUE, conf.method ="boot", nsim =1000)tbl_regression(mod_max_score_status_sex, intercept =TRUE, label =list(sex ~"Sex", migratory_status ~"Status"))# run rptR to obtain repeatabilities of random effectsrpt_mod_max_score_status_sex <-rpt(max_breeding_score ~ sex + migratory_status + (1| rings_comb),grname =c("rings_comb", "Fixed"),data = ind_breeding_scores_status,datatype ="Gaussian",nboot =1000, npermut =1000, ratio =TRUE,adjusted =TRUE, ncores =4, parallel =TRUE)# run partR2 on each model to obtain marginal R2, parameter estimates, and beta# weightsR2m_mod_max_score_status_sex <-partR2(mod_max_score_status_sex,partvars =c("sex", "migratory_status"),R2_type ="marginal",nboot =1000,CI =0.95,max_level =1)R2c_mod_max_score_status_sex <-partR2(mod_max_score_status_sex,partvars =c("sex", "migratory_status"),R2_type ="conditional",nboot =1000,CI =0.95,max_level =1)stats_mod_max_score_status_sex <-list(mod = mod_max_score_status_sex,tidy = tidy_mod_max_score_status_sex,rptR = rpt_mod_max_score_status_sex,partR2m = R2m_mod_max_score_status_sex,partR2c = R2c_mod_max_score_status_sex,data = ind_breeding_scores)# saveRDS(stats_mod_max_score_status_sex, file = "out/stats_mod_max_score_status_sex.rds")#### Table of effect sizes ----# Retrieve sample sizessample_sizes <- ind_breeding_scores_status %>%ungroup() %>%summarise(Year =n_distinct(season),Individual =n_distinct(rings_comb),Observations =nrow(.))sample_sizes <-as.data.frame(t(as.data.frame(sample_sizes))) %>%rownames_to_column("term") %>%rename(estimate = V1) %>%mutate(stat ="n")# clean model component namesmod_comp_names <-data.frame(comp_name =c("Male (Sex)","Resident (Migratory Status)","Total Marginal \U1D479\U00B2","Sex","Migratory Status","Total Conditional \U1D479\U00B2","Individual","Residual","Individual","Residual","Years","Individuals","Observations"))# Fixed effect sizes (non-standardized)fixefTable <- stats_mod_max_score_status_sex$tidy %>% dplyr::filter(effect =="fixed") %>% dplyr::select(term, estimate, conf.low, conf.high) %>%as.data.frame() %>%mutate(stat ="fixed")# Fixed effect sizes (standardized)fixef_bw_Table <- stats_mod_max_score_status_sex$partR2m$BW %>%as.data.frame() %>%mutate(stat ="fixed_bw") %>%rename(conf.low = CI_lower,conf.high = CI_upper)# Semi-partial R2 estimatesR2Table <-bind_rows(stats_mod_max_score_status_sex$partR2m$R2, stats_mod_max_score_status_sex$partR2c$R2[1,]) %>% dplyr::select(term, estimate, CI_lower, CI_upper) %>%as.data.frame() %>%mutate(stat ="partR2") %>%rename(conf.low = CI_lower,conf.high = CI_upper)# Random effects variancesranefTable <- stats_mod_max_score_status_sex$tidy %>% dplyr::filter(effect =="ran_pars") %>% dplyr::select(group, estimate, conf.low, conf.high) %>%as.data.frame() %>%mutate(stat ="rand") %>%rename(term = group) %>%mutate(estimate = estimate^2,conf.high = conf.high^2,conf.low = conf.low^2)# Adjusted repeatabilitiescoefRptTable <- stats_mod_max_score_status_sex$rptR$R_boot %>% dplyr::select(-Fixed) %>%mutate(residual =1-rowSums(.)) %>%apply(., 2, function(x) c(mean (x), quantile (x, prob =c(0.025, 0.975)))) %>%t() %>%as.data.frame() %>%rownames_to_column("term") %>%rename(estimate = V1,conf.low =`2.5%`,conf.high =`97.5%`) %>%mutate(stat ="RptR")# Store all parameters into a single table and clean it upallCoefs_mod <-bind_rows(fixef_bw_Table, R2Table, ranefTable, coefRptTable, sample_sizes) %>%bind_cols(., mod_comp_names) %>%mutate(coefString =ifelse(!is.na(conf.low),paste0("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"),NA),effect =c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)),rep("Partitioned \U1D479\U00B2", nrow(R2Table)),rep("Random effects \U1D70E\U00B2", nrow(ranefTable)),rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)),rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>% dplyr::select(effect, everything())``````{r}# draw gt tablemod_max_score_table_status_sex <- allCoefs_mod %>% dplyr::select(effect, comp_name, estimate, coefString) %>%gt(rowname_col ="row",groupname_col ="effect") %>%cols_label(comp_name =html("<i>Banded Dotterel rufous band score</i>"),estimate ="Mean estimate",coefString ="95% confidence interval") %>%fmt_number(columns =c(estimate),rows =1:10,decimals =2,use_seps =FALSE) %>%fmt_number(columns =c(estimate),rows =11:13,decimals =0,use_seps =FALSE) %>%sub_missing(columns =1:4,missing_text ="") %>%cols_align(align ="left",columns =c(comp_name)) %>%tab_options(row_group.font.weight ="bold",row_group.background.color =brewer.pal(9,"Greys")[3],table.font.size =12,data_row.padding =3,row_group.padding =4,summary_row.padding =2,column_labels.font.size =14,row_group.font.size =12,table.width =pct(80))mod_max_score_table_status_sexplot(allEffects(mod_max_score_status_sex))```### within-individual moult dynamics (sex)```{r}# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in sex while accounting for sex differences in maximum rufous band score.ind_prop_molt_scores <- dat %>%mutate(score =as.numeric(score)) %>%left_join(., select(ind_breeding_scores, -sex), by =c("rings_comb","season")) %>%filter(!is.na(max_breeding_score)) %>%mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%select(rings_comb, season, date_J, sex, score, max_breeding_score, prop_molt_score) %>%arrange(rings_comb, season, date_J) %>%distinct() %>%filter(!is.na(prop_molt_score))# Assess sample sizes of each sexind_prop_molt_scores %>%group_by(sex) %>%summarise(n_distinct(rings_comb))# Step 1 and 2: Standardize max_breeding_score within each sex groupind_prop_molt_scores <- ind_prop_molt_scores %>%group_by(sex) %>%mutate(std_max_breeding_score = (max_breeding_score -mean(max_breeding_score, na.rm =TRUE)) /sd(max_breeding_score, na.rm =TRUE)) %>%ungroup()# mixed effects binomial model comparing sex and date effect on the changes in moult scoresmod1 <- lme4::glmer(prop_molt_score ~ date_J * sex + std_max_breeding_score +# date_J * sex + (1| rings_comb) + (1|season),data = ind_prop_molt_scores, family = binomial,control =glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =1e+06)))plot(allEffects(mod1))# strong date effect, but no sex effecttbl_regression(mod1, intercept =TRUE, label =list(date_J ~"Date", sex ~"Sex"))# extract predicted trendspr <- ggeffects::predict_response(mod1, c("date_J [30:293]", "sex"))# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian datedates_for_plot <-data.frame(date =as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),date_J =c(1:365))# join the back-transformed dates to model fitsmod1_fits <-as.data.frame(pr) %>%rename(date_J = x,sex = group) %>%left_join(., dates_for_plot, by ="date_J")# plot the modelggplot() +geom_line(data = mod1_fits, aes(x = date, y = predicted, color = sex)) +geom_ribbon(data = mod1_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = sex),lwd =1, alpha =0.25) + luke_theme +theme(legend.position =c(0.3, 0.2),legend.justification =c(1, 0),strip.background =element_blank(),axis.title.x =element_blank(),axis.text.x =element_text(size =10, angle =45, hjust =1, vjust =1)) +ylab("proportion intact of individual-based maximum rufous band") +scale_x_date(date_labels ="%B", expand =c(0.01, 0.01), date_breaks ="1 month") +scale_colour_brewer(palette ="Dark2", direction =-1,name ="Sex",labels =c("Female (N = 49)", "Male (N = 35)")) +scale_fill_brewer(palette ="Dark2", direction =-1,name ="Sex",labels =c("Female (N = 49)", "Male (N = 35)"))```### within-individual moult dynamics (migratory status)```{r}# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in migratory status.ind_prop_molt_scores <- dat %>%mutate(score =as.numeric(score)) %>%left_join(., select(ind_breeding_scores, -sex), by =c("rings_comb","season")) %>%filter(!is.na(max_breeding_score)) %>%mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%select(rings_comb, season, date_J, sex, migratory_status, score, max_breeding_score, prop_molt_score) %>%arrange(rings_comb, season, date_J) %>%distinct() %>%filter(!is.na(prop_molt_score)) %>%filter(migratory_status %in%c("R", "M"))# Assess sample sizes of each sexind_prop_molt_scores %>%group_by(sex, migratory_status) %>%summarise(n_distinct(rings_comb))# mixed effects binomial model comparing sex and date effect on the changes in moult scoresmod1_status <- lme4::glmer(prop_molt_score ~ date_J * migratory_status +# date_J * sex + (1| rings_comb) + (1| season),data = ind_prop_molt_scores, family = binomial,control =glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =1e+06)))# strong date effect, but no status effecttbl_regression(mod1_status, intercept =TRUE, label =list(date_J ~"Date", migratory_status ~"Migratory status"))# extract predicted trendspr <- ggeffects::predict_response(mod1_status, c("date_J [30:293]", "migratory_status"))# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian datedates_for_plot <-data.frame(date =as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),date_J =c(1:365))# join the back-transformed dates to model fitsmod1_fits <-as.data.frame(pr) %>%rename(date_J = x,migratory_status = group) %>%left_join(., dates_for_plot, by ="date_J")# plot the modelggplot() +geom_line(data = mod1_fits, aes(x = date, y = predicted, color = migratory_status)) +geom_ribbon(data = mod1_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = migratory_status),lwd =1, alpha =0.25) + luke_theme +theme(legend.position =c(0.3, 0.2),legend.justification =c(1, 0),strip.background =element_blank(),axis.title.x =element_blank(),axis.text.x =element_text(size =10, angle =45, hjust =1, vjust =1)) +ylab("proportion intact of individual-based maximum rufous band") +scale_x_date(date_labels ="%B", expand =c(0.01, 0.01), date_breaks ="1 month") +scale_colour_brewer(palette ="Set1", direction =-1,name ="Migratory status",labels =c("Resident (N = 16)", "Migratory (N = 5)")) +scale_fill_brewer(palette ="Set1", direction =-1,name ="Migratory status",labels =c("Resident (N = 16)", "Migratory (N = 5)"))```### within-individual moult dynamics (age)```{r}# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in migratory status.ind_prop_molt_scores <- dat %>%mutate(score =as.numeric(score)) %>%left_join(., select(ind_breeding_scores, -sex), by =c("rings_comb","season")) %>%filter(!is.na(max_breeding_score)) %>%mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%select(rings_comb, season, date_J, sex, age, score, max_breeding_score, prop_molt_score) %>%arrange(rings_comb, season, date_J) %>%distinct() %>%filter(!is.na(prop_molt_score)) %>%filter(!is.na(age))# Assess sample sizes of each sexind_prop_molt_scores %>%group_by(age) %>%summarise(n_distinct(rings_comb))# mixed effects binomial model comparing sex and date effect on the changes in moult scoresmod1_age <- lme4::glmer(prop_molt_score ~ date_J + age +# date_J * sex + (1| rings_comb),# + (1 | season),data = ind_prop_molt_scores, family = binomial)# control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))# strong date effect, but no status effecttbl_regression(mod1_age, intercept =TRUE, label =list(date_J ~"Date", age ~"Age"))# extract predicted trendspr <- ggeffects::predict_response(mod1_age, c("date_J [30:293]", "age"))# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian datedates_for_plot <-data.frame(date =as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),date_J =c(1:365))# join the back-transformed dates to model fitsmod1_fits <-as.data.frame(pr) %>%rename(date_J = x,age = group) %>%left_join(., dates_for_plot, by ="date_J")# plot the modelggplot() +geom_line(data = mod1_fits, aes(x = date, y = predicted, color = age)) +geom_ribbon(data = mod1_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = age),lwd =1, alpha =0.25) + luke_theme +theme(legend.position =c(0.3, 0.2),legend.justification =c(1, 0),strip.background =element_blank(),axis.title.x =element_blank(),axis.text.x =element_text(size =10, angle =45, hjust =1, vjust =1)) +ylab("proportion intact of individual-based maximum rufous band") +scale_x_date(date_labels ="%B", expand =c(0.01, 0.01), date_breaks ="1 month", limits =c(as.Date("2021-10-01"), as.Date("2022-05-01"))) +scale_colour_brewer(palette ="Spectral", direction =-1,name ="Age",labels =c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", "4 (N = 13", "5 (N = 4)", "6 (N = 2)")) +scale_fill_brewer(palette ="Spectral", direction =-1,name ="Age",labels =c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", "4 (N = 13", "5 (N = 4)", "6 (N = 2)"))```### 13-Nov-24: meeting notes (Bart, Luke, and Bashar via zoom)- add the breeding data to the analysis: for each individual include the date of last breeding (i.e., last date seen with nest or brood) as a predictor of molt timing. Also do an analysis relating the timing of last breeding to the age of the individual. Bashar will send Luke the breeding data he received from Ted.- drop pre-October data from the molt timing analysis.- next steps: Bashar to eventually write up this study as a report (which could eventually be used as his MSc thesis if he wants). Priority will be on him writing up a proposal for his thesis (which has a concrete deadline and is graded). Bashar will present the study at our weekly Monday meeting in the first few months of 2025.